Построение модели и прогнозирование на июль 2018
In [1]:
# basic stuff
import pandas as pd
import numpy as np
from itertools import product
import holidays
import pickle
pd.set_option('display.max_columns', 500)
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
# time stuff
from datetime import datetime
hour = pd.DateOffset(hours=1)
# modeL stuff
from sklearn.metrics import pairwise_distances_argmin_min, r2_score
from sklearn.model_selection import TimeSeriesSplit
import xgboost as xgb
from xgboost import plot_importance, XGBRegressor
# geo stuff
import geopandas as gpd
import json
# plotting stuff
import plotly
import plotly.graph_objs as go
from plotly.subplots import make_subplots
import plotly.io as pio
from plotly.offline import iplot
import matplotlib.pyplot as plt
In [2]:
pio.templates['rockwell'] = go.layout.Template(
layout=go.Layout(
font=dict(
family='Rockwell',
size=18,
color='#2A5674'
),
title={'yanchor': 'middle'},
)
)
pio.templates.default = "plotly_white+rockwell"
In [3]:
pink= '#e2d9e2'
blue = '#6785be'
red = '#8e2c50'
orange = '#ba6657'
beige = '#ceac94'
purple = '#421257'
In [4]:
zones_gdf = gpd.read_file('NYC_TAXI_data/other_data/taxi_zones.geojson').loc[:,['zone', 'OBJECTID', 'borough','geometry']]
zones_gdf.columns = ['zone_name', 'zone_id', 'borough','geometry']
zones_gdf['zone_id'] = zones_gdf['zone_id'].astype(str)
NY_center_lat = (40.49612+40.91553)/2
NY_center_lon = (-74.25559-73.70001)/2
In [5]:
zones_gdf.sample(5)
Out[5]:
In [6]:
zones_geojson = json.loads(zones_gdf.to_json())
layout = go.Layout(mapbox_style='carto-positron',
mapbox_zoom=9.5,
mapbox_center = {'lat': NY_center_lat, 'lon': NY_center_lon},
hoverlabel=dict(
bgcolor="white",
font_size=12,
font_family='Rockwell'),
margin={"r":0,"t":50,"l":0,"b":0},
title='New York taxi zones',
height=600, width=800
)
data = go.Choroplethmapbox(geojson=zones_geojson,
locations=zones_gdf['zone_id'],
z=zones_gdf['zone_id'],
hovertext=zones_gdf['zone_name'],
hovertemplate='<b>Zone name</b>: <b>%{hovertext}</b>'+
'<br><b>Zone ID </b>: %{location}'+
"<extra></extra>",
showlegend=False,
autocolorscale=False,
colorscale='twilight',
showscale=False,
marker_opacity=0.8,
marker_line_width=0.1
)
taxi_zones = go.Figure(data=data, layout=layout)
In [7]:
taxi_zones.write_html("exp/interactive_taxi_zones.html")
In [8]:
taxi_zones.show()
# при наведении на зону можно увидеть ее название и ID
2. Временные ряды¶
Для того, чтобы решить поставленную задачу, сырые данные необходимо агрегировать по часам и районам. Агрегированные данные будут представляют собой почасовые временные ряды с количеством поездок из каждой зоны:
In [9]:
aggregated_data = pd.read_csv('NYC_TAXI_aggregated_data/pu_agg_data.csv',
index_col=['pickup_datetime'], parse_dates=['pickup_datetime'])
aggregated_data.index.name = 'pickup_datetime'
aggregated_data.columns.name = 'zone_id'
In [10]:
aggregated_data.head()
Out[10]:
Нарисуем несколько временных рядов, чтобы посмотреть, как они вообще выглядят:
In [11]:
timerows_to_plot = ['7', '161', '261']
timerows_colors = {timerows_to_plot[0]: purple,
timerows_to_plot[1]: blue,
timerows_to_plot[2]: orange
}
interactive_some_timerows = go.Figure()
for zone in timerows_to_plot:
df = aggregated_data[[zone]].reset_index()
x = df['pickup_datetime']
y = df[str(zone)]
interactive_some_timerows.add_trace(go.Scatter(x=x, y=y,
mode='lines',
name=str(zone),
line=dict(width=2,
color=timerows_colors[zone])))
interactive_some_timerows = interactive_some_timerows.update_xaxes(rangeslider_visible=True)
interactive_some_timerows = interactive_some_timerows.update_layout(title='Timerows for selected taxi zones',
autosize=True)
In [12]:
interactive_some_timerows.write_html("exp/interactive_some_timerows.html")
In [13]:
interactive_some_timerows.show()
# можно выбрать интересующий период на шкале снизу, а сверху все будет в увеличенном виде!
Сильно выражена суточная сезонность, а так же недельная и годовая. Ряды разные.
Обучение и тест
Нужно обозначить дату, разделяющую обучение и тест. Обучаться будем до 2018/06/30 23:00, а предсказывать будем на июль 2018 года.
In [14]:
split_date = datetime(2018,6,30,23) # время, разделяющее train и test
2.1 Отбор рядов и визуализация спроса на такси на карте¶
Не все ряды одинаково полезны, то есть не все зоны одинаково востребованы. Прогнозировать спрос имеет смысл в популярных зоных, чтобы понимать, сколько в зоне должно быть доступных машин в каждый час.
Посчитаем среднее число поездок из каждой зоны в час и раскрасим регионы в соответствии с полученными значениями:¶
In [15]:
mean_counts = pd.DataFrame(aggregated_data.loc[:split_date].mean(), columns=['mean_count']).reset_index()
zones_gdf_with_means = zones_gdf.merge(mean_counts, how='left', on='zone_id')#.fillna(0)
In [16]:
zones_gdf_with_means.head()
Out[16]:
In [17]:
zones_geojson = json.loads(zones_gdf_with_means.to_json())
layout = go.Layout(
mapbox_style='carto-positron',
mapbox_zoom=9,
mapbox_center = {'lat': NY_center_lat, 'lon': NY_center_lon},
hoverlabel=dict(
bgcolor="white",
font_size=12,
font_family='Rockwell'),
margin={"r":0,"t":50,"l":0,"b":0},
title='New York taxi zones colored by mean trip count',
height=600, width=600
)
data = go.Choroplethmapbox(geojson=zones_geojson,
locations=zones_gdf_with_means['zone_id'],
z=zones_gdf_with_means['mean_count'],
hovertext=zones_gdf_with_means['zone_name'],
hovertemplate='<b>Zone name</b>: <b>%{hovertext}</b>'+
'<br><b>Zone ID </b>: %{location}'+
'<br><b>Mean trip count </b>: %{z}<br>'+
"<extra></extra>",
showlegend=False,
autocolorscale=False,
colorscale='BuPu',
marker_opacity=0.8, marker_line_width=0.1
)
mean_colored_zones = go.Figure(data=data, layout=layout)
In [18]:
mean_colored_zones.write_html("exp/mean_colored_zones.html")
In [19]:
mean_colored_zones.show()
Работать будем с зонами, из которых среднее число поездок в час не меньше 5, таких зон получется 83:
In [20]:
popular_zones = list(aggregated_data.loc[:, aggregated_data.mean() >= 5].columns)
print(len(popular_zones))
In [21]:
R = len(popular_zones)
In [22]:
agg_data_pop = aggregated_data[popular_zones].copy()
agg_data_pop.to_csv('NYC_TAXI_aggregated_data/pu_agg_data_pop.csv', index=True)
agg_data_pop.head()
Out[22]:
In [23]:
del aggregated_data
In [24]:
popmean_counts = pd.DataFrame(agg_data_pop.mean(), columns=['mean_count']).reset_index()
popzones_with_means = zones_gdf.merge(popmean_counts, how='left', on='zone_id')
popzones_with_means.head()
Out[24]:
In [25]:
popzones_geojson = json.loads(popzones_with_means.to_json())
layout = go.Layout(mapbox_style='carto-positron',
mapbox_zoom=9,
mapbox_center = {'lat': NY_center_lat, 'lon': NY_center_lon},
hoverlabel=dict(
bgcolor="white",
font_size=12,
font_family='Rockwell'),
margin={"r":0,"t":50,"l":0,"b":0},
title='Most popular New York taxi zones<br>colored by mean trip count',
height=600, width=600)
data = go.Choroplethmapbox(geojson=popzones_geojson,
locations=popzones_with_means['zone_id'],
z=popzones_with_means['mean_count'],
hovertext=popzones_with_means['zone_name'],
hovertemplate='<b>Zone name</b>: <b>%{hovertext}</b>'+
'<br><b>Zone ID </b>: %{location}'+
'<br><b>Mean trip count </b>: %{z}<br>'+
"<extra></extra>",
showlegend=False,
autocolorscale=False,
colorscale='BuPu',
marker_opacity=0.8, marker_line_width=0.1
)
mean_colored_pop_zones = go.Figure(data=data, layout=layout)
In [26]:
mean_colored_pop_zones.show()
2.2. Кластеризация рядов¶
Так как подбирать параметры для отдельных моделей каждой географической зоны относительно долго, то кластеризуем ряды и подбирать параметры модели будем уже для каждого кластера отдельно.
Кластеризацию делаю методом KMeans.
In [27]:
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.decomposition import PCA
In [28]:
# стандартизуем временные ряды
scaler = StandardScaler()
data_for_clusterization = pd.DataFrame(scaler.fit_transform(agg_data_pop.loc[:split_date]),
columns=agg_data_pop.loc[:split_date].columns,
index=agg_data_pop.loc[:split_date].index).T
In [29]:
clusterization = make_subplots(rows=2, cols=3, subplot_titles=('n_clusters = 2',
'n_clusters = 3',
'n_clusters = 4',
'n_clusters = 5',
'n_clusters = 6',
'n_clusters = 7',
))
for n_clusters, (row, col) in zip(range(2,8), product([1,2],[1,2,3])):
# кластеризуем
clusterer = KMeans(n_clusters=n_clusters, random_state=21)
labels = clusterer.fit_predict(data_for_clusterization)
centers = clusterer.cluster_centers_
# уменьшим размерность
reductor = PCA(2, random_state=21)
dots = reductor.fit_transform(data_for_clusterization)
dots = pd.DataFrame(dots, columns=['x','y'])
dots['labels'] = labels.astype(float)
# добавим 2D-график
clusterization.add_trace(go.Scatter(x=dots['x'], y=dots['y'], mode='markers',
marker=dict(
size=10,
line_width=0.5,
color=dots['labels'],
colorscale='Viridis',
showscale=False,
),
text=dots['labels'],
hovertemplate='<b>Cluster #%{text}</b>'+'<extra></extra>'
),
row=row, col=col)
clusterization.update_xaxes(title_text='silhouette_score {:.2f}'.format(silhouette_score(data_for_clusterization, labels)),
title_font=dict(size=14, family='Rockwell'),
row=row, col=col)
clusterization = clusterization.update_layout(height=600, width=800, title_text='Clusterization of popular zones', showlegend=False)
clusterization = clusterization.update_xaxes(showticklabels=False)
clusterization = clusterization.update_yaxes(showticklabels=False)
In [30]:
clusterization.write_html("exp/clusterization.html")
In [31]:
clusterization.show()
# кластеризация стандартизованных временных рядов методом KMeans, подсчет коэффициента силуэта
Делить будем на 3 кластера.
In [32]:
n_clusters = 3
clusterer = KMeans(n_clusters=n_clusters, random_state=21)
labels = clusterer.fit_predict(data_for_clusterization)
centers = clusterer.cluster_centers_
In [33]:
zone_cluster_labels = {zone:label for zone, label in zip(popular_zones, labels)}
Найдем ряды, ближайшие к центрам кластеров:
In [34]:
# центры кластеров
closest, _ = pairwise_distances_argmin_min(centers, data_for_clusterization)
centroid_timerows = agg_data_pop.iloc[:,closest].copy()
centroid_names = centroid_timerows.columns.values # типичные ряды для каждого кластера
with open('XGBoost/centroid_names.pickle','wb') as file:
pickle.dump(centroid_names, file)
In [35]:
del data_for_clusterization
In [36]:
print('Ближайшие к центрам кластеров ряды: {}'.format(centroid_names))
In [37]:
timerows_colors = {centroid_names[0]: purple,
centroid_names[1]: blue,
centroid_names[2]: orange,
}
centroid_timerows = go.Figure()
for zone in centroid_names:
x = agg_data_pop[[zone]].index
y = agg_data_pop[[zone]][str(zone)]
centroid_timerows.add_trace(go.Scatter(x=x, y=y, mode='lines',
name=str(zone),
line=dict(width=2, color=timerows_colors[zone])))
centroid_timerows = centroid_timerows.update_xaxes(rangeslider_visible=True)
centroid_timerows = centroid_timerows.update_layout(title='Timerows for "centroid" taxi zones',
autosize=True)
In [38]:
centroid_timerows.write_html("exp/centroid_timerows.html")
In [39]:
centroid_timerows.show()
3. Подготовка выборок для разных задач прогнозирования¶
Для каждой из шести задач прогнозирования $\hat{y}_{T+i|T}, i=1,\dots,6$ будем формировать выборки, в которых отклик = $y_{T+i}$ при всевозможных значениях $T$ (то есть двигаем значения ряда "назад").
In [40]:
shifted_data = {}
for shift in range(1,7):
shifted_data[shift] = agg_data_pop.shift(-shift).dropna()
4. Построение признаковых описаний рядов¶
Теперь к временному ряду каждого региона добавим признаки. Признаки берутся из:
- Значений времени ряда:
- год, месяц, день месяца, день недели, час
- праздник/обычный день
- сезонные гармоники
Признаки, зависящие от значений ряда (авторегрессионные):
- cреднее количество поездок за предыдущие 12/24/48/168 часов
- суммарное количество поездок за предыдущие 12/24/48/168 часов
- среднее число поездок в такой час/день/день недели/выходной
- количество поездок из рассматриваемого района в моменты времени $y_{T-1}, ... , y_{T-K}$
- количество поездок из рассматриваемого района в моменты времени $y_{T-24}, y_{T-48}, ..., y_{T-24*K_d}$
Статистической информации о регионе из сырых данных (эти признаки напрямую не зависят от значений целевой переменной):средняя длительность поездоксреднее количество пассажировсреднее расстояние по счётчикудоли географических зон, в которые совершаются поездкидоли поездок, совершаемых по тарифам каждого из типовдоли способов оплаты поездоксредняя стоимость поездокдоли провайдеров данных
4.1. Метаданные для создания признаков¶
Загружаем размеченные по регионам записи о поездках - marked_data за каждый месяц. Я не могу хранить такие данные в памяти, поэтому надо подумать здесь о БД!
4.2. Функция, для добавления признаков в датафрейм с временным рядом¶
Признаков делаем много, потом проведем отбор.
In [41]:
def add_regressors(timeseries, K_s, K_h, K_d, meta_data=None):
timerow = timeseries.copy() # сюда будем записывать регрессоры
timerow['datetime'] = timerow.index # вспомогательный столбец
### 1.1 Признаки, зависящие от datetime
# год, месяц, день месяца, день недели, час
timerow['year'] = timerow['datetime'].dt.year
timerow['month'] = timerow['datetime'].dt.month
timerow['quarter'] = timerow['datetime'].dt.quarter
timerow['day'] = timerow['datetime'].dt.day
timerow['hour'] = timerow['datetime'].dt.hour
timerow['weekday'] = timerow['datetime'].dt.dayofweek
timerow['dayofyear'] = timerow['datetime'].dt.dayofyear
timerow['weekofyear'] = timerow['datetime'].dt.week
timerow['weekend'] = (timerow['datetime'].dt.weekday>=5).astype(int)
# праздник/обычный день
us_holidays = holidays.UnitedStates() # данные об американских праздниках
timerow['holiday'] = np.array([date in us_holidays for date in timerow['datetime']]).astype(int)
# сезонные гармоники
T=len(timerow.index)
for i in range(1,K_s+1):
timerow['dsin_'+str(i)] = np.sin(np.linspace(1,T,T)*((2*np.pi*i)/24.))
timerow['dcos_'+str(i)] = np.cos(np.linspace(1,T,T)*((2*np.pi*i)/24.))
timerow['wsin_'+str(i)] = np.sin(np.linspace(1,T,T)*((2*np.pi*i)/168.))
timerow['wcos_'+str(i)] = np.cos(np.linspace(1,T,T)*((2*np.pi*i)/168.))
timerow['ysin_'+str(i)] = np.sin(np.linspace(1,T,T)*((2*np.pi*i)/8766.))
timerow['ycos_'+str(i)] = np.cos(np.linspace(1,T,T)*((2*np.pi*i)/8766.))
### 1.2 Признаки, зависящие от предыдущих значений ряда
#cреднее количество поездок за предыдущие 12/24/48/168 часов
timerow['prev_mean_12'] = timerow['real_counts'].rolling(window=12).mean()
timerow['prev_mean_24'] = timerow['real_counts'].rolling(window=24).mean()
timerow['prev_mean_48'] = timerow['real_counts'].rolling(window=48).mean()
timerow['prev_mean_168'] = timerow['real_counts'].rolling(window=168).mean()
# суммарное количество поездок за предыдущие 12/24/48/168 часов
timerow['prev_sum_12'] = timerow['real_counts'].rolling(window=12).sum()
timerow['prev_sum_24'] = timerow['real_counts'].rolling(window=24).sum()
timerow['prev_sum_48'] = timerow['real_counts'].rolling(window=48).sum()
timerow['prev_sum_168'] = timerow['real_counts'].rolling(window=168).sum()
# среднее число поездок в такой час/день/день недели/выходной
timerow['hour_mean'] = timerow.groupby(by=['hour'])['real_counts'].transform('mean')
timerow['day_mean'] = timerow.groupby(by=['day'])['real_counts'].transform('mean')
timerow['dayofweek_mean'] = timerow.groupby(by=['weekday'])['real_counts'].transform('mean')
timerow['weekend_mean'] = timerow.groupby(by=['weekend'])['real_counts'].transform('mean')
# количество поездок из рассматриваемого района в моменты времени $y_{T-1}, ... , y_{T-K}$
for i in range(1,K_h+1):
timerow['prev_count_'+str(i)] = timerow['real_counts'].shift(i)
# количество поездок из рассматриваемого района в моменты времени $y_{T-24}, y_{T-48}, ..., y_{T-24*K_d}$
for i in range(1,K_d+1):
timerow['prev_count_'+str(i*24)] = timerow['real_counts'].shift(i*24)
### 1.3 Признаки, зависящие от метаданных географической зоны
# идентификатор зоны
zone = timeseries.columns.name
timerow['zone'] = int(zone)
#a = meta_data.groupby(by=['pickup_datetime']).agg(np.mean)
#timerow['mean_passenger_count'] = a['passenger_count']# среднее количество пассажиров
#timerow['mean_trip_distance'] = a['trip_distance']# среднее расстояние по счётчику
#timerow['mean_trip_duration'] = a['trip_duration']# средняя длительность поездок
#timerow['mean_total_cost'] = a['total_cost']# средняя стоимость поездок
# доли способов оплаты поездок
#for x in np.unique(meta_data['payment_type']):
# timerow['payment_type_{}_ratio'.format(int(x))] = meta_data[meta_data['payment_type'] == x]['payment_type'].groupby(by=['pickup_datetime']).count()\
# / timerow['real_counts']
# доли провайдеров данных
#for x in np.unique(meta_data['vendorID']):
# timerow['vendor_{}_ratio'.format(int(x))] = meta_data[meta_data['vendorID'] == x]['vendorID'].groupby(by=['pickup_datetime']).count()\
# / timerow['real_counts']
timerow.drop('datetime', axis=1, inplace=True)
timerow.fillna(0, inplace=True)
return timerow
In [42]:
#по умолчанию
K_s, K_h, K_d = (30,14,7) #сделаем много авторегрессионных признаков, а потом устроим отбор
In [43]:
zone = centroid_names[0]
timeseries = agg_data_pop.loc[:, [zone]]
timeseries.columns = ['real_counts']
timeseries.columns.name = zone
#zone_data = marked_data[marked_data['PU_zone_id'].astype(str) == zone]
timeseries_with_features = add_regressors(timeseries, K_s, K_h, K_d)
timeseries_with_features.columns.name = zone
После применения функции добавления признаков к ряду получается такой датафрейм:
In [44]:
timeseries_with_features.sample(5)
Out[44]:
Данные подготовили, дальше - модель. Я выбрала XGBoost (быстро, удобно, качественно).
5. Отбор признаков для XGBoost¶
Некоторые признаки могут быть несущественными или незначимыми для целевой переменной, и их включение в модель приводит к увеличению сложности модели, которую сложнее интерпретировать, а так же больше количество признаков увеличивает время обучения модели и уменьшает надежность предсказаний.
In [45]:
# define custom class to fix bug in xgboost 1.1.0
class myXGBR(XGBRegressor):
@property
def coef_(self):
return None
Отбор признаков проводила следующим образом. Для центроида каждого кластера по кросс-валидации для временных рядов [1] оценивала качество модели, обученной на всех созданных признаках, по метрике r2_score, и фильтровала признаки по значению features_importances (XGBR.featureimportances) с порогом 0.01. Затем обучала модель на отобранных признаках и снова считала метрику r2. Конечно, это все отдельно для каждой прогнозирующей задачи.
Посмотрим, как признаки используются в моделях для каждого кластера. F-score на графиках ниже - частота появления признака в деревьях. Некоторые признаки почти не используются. Почистим от них. Threshold возьмем равным 0.01.
In [46]:
def cv_r2_score(n_splits, X_train, y_train, params={'objective': 'reg:squarederror', 'subsample': 1}):
tscv = TimeSeriesSplit(n_splits=3)
scores = []
importances = []
for tr_ix, cv_ix in tscv.split(X_train):
XGBR = myXGBR()
XGBR.fit(X_train.iloc[tr_ix], y_train.iloc[tr_ix])
y_pred = XGBR.predict(X_train.iloc[cv_ix])
y_test = y_train.iloc[cv_ix].values
scores.append(r2_score(y_test, y_pred))
importances.append(XGBR.feature_importances_)
score = np.mean(scores) # функция возвращает средне по cv значение r2
feature_importances = np.array(importances).mean(axis=0) # и оценку важности признаков
return score, feature_importances
In [47]:
%%time
features_dict = {}
dict_of_r2 = {}
for shift in range(1,7):
data = shifted_data[shift]
print('Модель T+{}; {}'.format(shift, datetime.now().strftime('%H:%M:%S')))
features_dict[shift] = {}
dict_of_r2[shift] = {}
fig, ax = plt.subplots(1, 3, figsize=(15,5), dpi=500)
for cluster_label, zone in enumerate(centroid_names):
print('Кластер {}, зона {}; {}'.format(cluster_label, zone, datetime.now().strftime('%H:%M:%S')))
# временной ряд региона
timeseries = data.loc[:, [zone]]
timeseries.columns = ['real_counts']
timeseries.columns.name = zone
# метаданные по региону
#zone_data = marked_data[marked_data['PU_zone_id'] == zone].shift(-shift).dropna()
# добавляем признаки
timeseries_with_features = add_regressors(timeseries, K_s, K_h, K_d)
timeseries_with_features.columns.name = zone
# отделяем train
timeseries_train = timeseries_with_features.loc[:split_date].copy()
X_train = timeseries_train.drop(['real_counts'], axis=1)
y_train = timeseries_train[['real_counts']]
# оцениваем по кросс-валидации модель, обученную на всех признаках с параметрами по умолчанию
r2_, feature_importances = cv_r2_score(4, X_train, y_train)
dict_of_r2[shift][cluster_label] = [r2_]
print(u'r2_score для всех признаков: {:.2f}'.format(r2_))
# отбираем признаки из модели
features_df = pd.DataFrame({'features': X_train.columns,
'importances': feature_importances}).sort_values('importances',
ascending=False)
selected_features = features_df[features_df['importances'] > 0.01]['features'].values
features_dict[shift][cluster_label] = selected_features
print('{} наиболее релевантных признаков'.format(len(selected_features)))
# оцениваем по кросс-валидации модель, обученную на выбранных признаках с параметрами по умолчанию
r2_, _ = cv_r2_score(4, X_train[selected_features], y_train)
dict_of_r2[shift][cluster_label].append(r2_)
print(u'r2_score для отобранных признаков: {:.2f}'.format(r2_))
# визуализизуем релевантность отобранных признаков
XGBR = myXGBR()
XGBR.fit(X_train[selected_features], y_train)
ax_ = ax[cluster_label]
ax_ = plot_importance(XGBR, ax=ax_)
ax_.set_title(u'Релевантность признаков для кластера {}'.format(cluster_label), size=12)
ax_.set_ylabel('')
print('===')
plt.subplots_adjust(wspace = .6, hspace = .3)
plt.suptitle(u'Отбор признаков для модели T+{}'.format(shift))
plt.show()
Даже если метрика не улучшилась, то уменьшение числа признаков ускоряет обучение.
In [48]:
with open('XGBoost/selected_features.pickle','wb') as file:
pickle.dump(features_dict, file)
In [49]:
# посмотрим, как изменилась метрика r2 после отбора признаков
shift = 3
r2_df = pd.DataFrame(dict_of_r2[shift]).T
r2_df.columns=['all_features/default_params',
'selected_features/default_params'
]
r2_df.index.name = 'cluster'
Например, так изменилась метрика R2 модели Т+3 при отборе признаков:
In [50]:
r2_df
# изменение метрики R2 модели Т+3 при отборе признаков
Out[50]:
6. Подбор гиперпараметров для XGBoost¶
Теперь надо подобрать гиперпараметры моделей каждого кластера. Будем оценивать модель с разными параметрами по метрике R2 на кросс-валидации через Bayesian Optimization Process (библиотека)
https://krasserm.github.io/2018/03/21/bayesian-optimization/ - разобраться при случае с этой статьей
In [51]:
from bayes_opt import BayesianOptimization
In [52]:
with open('XGBoost/selected_features.pickle', 'rb') as file:
features_dict = pickle.load(file)
In [53]:
# функция, которую будем оптимизировать
def cv_score_opt(learning_rate, max_depth, gamma, n_estimators, reg_alpha, reg_lambda):
xgb_params = {
'objective': 'reg:squarederror',
'subsample': 1,
'learning_rate': learning_rate,
'max_depth': int(max_depth),
'gamma': gamma,
'n_estimators': int(n_estimators),
'reg_alpha': reg_alpha,
'reg_lambda': reg_lambda
}
# Cross-validation для временных рядов
tscv = TimeSeriesSplit(n_splits=3)
scores = []
for tr_ix, cv_ix in tscv.split(X_train):
XGBR = myXGBR(**xgb_params)
XGBR.fit(X_train.iloc[tr_ix], y_train.iloc[tr_ix])
y_pred = XGBR.predict(X_train.iloc[cv_ix])
y_test = y_train.iloc[cv_ix].values
scores.append(r2_score(y_test, y_pred))
score = np.mean(scores)
return score
In [54]:
# сетка параметров
params_grid = {
'learning_rate': (0,1),
'max_depth': (1,100),
'gamma': (0,10),
'n_estimators': (10,300),
'reg_alpha': (0,1),
'reg_lambda': (0,1)
}
In [55]:
%%time
best_hyperparams = {}
for shift in range(1,7):
data = shifted_data[shift]
print(f'Модель T+{shift}; {datetime.now().strftime("%H:%M:%S")}')
best_hyperparams[shift] = {}
for cluster_label, zone in enumerate(centroid_names):
print(f'Кластер {cluster_label}, зона {zone}; {datetime.now().strftime("%H:%M:%S")}')
# временной ряд региона
timeseries = data.loc[:, [zone]]
timeseries.columns = ['real_counts']
timeseries.columns.name = zone
# метаданные по региону - не используем
#zone_data = marked_data[marked_data['PU_zone_id'] == zone].shift(-shift).dropna()
# добавляем признаки
timeseries_with_features = add_regressors(timeseries, K_s, K_h, K_d)
timeseries_with_features.columns.name = zone
# отделяем train и выбирыем нужные признаки
selected_features = features_dict[shift][cluster_label]
timeseries_train = timeseries_with_features.loc[:split_date].copy()
X_train = timeseries_train.drop(['real_counts'], axis=1)[selected_features]
y_train = timeseries_train[['real_counts']]
# инициализируем BO оптимизатор
optimizer = BayesianOptimization(cv_score_opt, params_grid, verbose=1)
# максимизируем r2_score
optimizer.maximize(n_iter=25, init_points=5, acq='ei')
# подобранные лучшие параметры
best_p = optimizer.max['params']
best_p['max_depth'] = np.ceil(best_p['max_depth'])
best_p['n_estimators'] = np.ceil(best_p['n_estimators'])
# запишем в словарь
xgb_params = {'objective': 'reg:squarederror',
'subsample': 1} # это константные параметры
xgb_params.update(best_p)
best_hyperparams[shift][cluster_label] = xgb_params
print(xgb_params)
# оцениваем по кросс-валидации модель, обученную на выбранных признаках с подобранными параметрами
r2_, _ = cv_r2_score(4, X_train, y_train, params=xgb_params)
print(u'r2_score для отобранных признаков c настроенными параметрами: {:.2f}'.format(r2_))
dict_of_r2[shift][cluster_label].append(r2_)
print('===')
In [56]:
with open('XGBoost/best_hyperparams.pickle','wb') as file:
pickle.dump(best_hyperparams, file)
In [57]:
# посмотрим, как изменилась метрика r2 после настройки гиперпараметров модели
shift = 3
r2_df = pd.DataFrame(dict_of_r2[shift]).T
r2_df.columns=['all_features/default_params',
'selected_features/default_params',
'selected_features/optimized_params'
]
r2_df.index.name = 'cluster'
In [58]:
r2_df
Out[58]:
7. Прогнозирование на июль 2018 года¶
Выбранными моделями построим для каждой географической зоны и каждого конца истории от 2018.06.30 23:00 до 2018.07.31 17:00 прогнозы на 6 часов вперёд; посчитаем ноутбуке ошибку прогноза по следующему функционалу:
$$Q_{july2018}=\frac{1}{R∗739∗6}\sum_{r=1}^{R}\sum_{T=2018.06.30 23:00}^{2018.07.31 17:00}\sum_{i=1}^{6}∣\widehat{y}^{r}_{T|T+i}−{y}^{r}_{T+i}| $$, где R - это число рядов (83 у нас), 739 - количество часов в прогнозируемом периоде.
In [59]:
# прогнозируем на июнь
july = pd.date_range(datetime(2018,6,30,23), datetime(2018,7,31,23), freq='H')
hour = pd.DateOffset(hours=1)
In [60]:
len(july)-6
Out[60]:
In [61]:
with open('XGBoost/selected_features.pickle', 'rb') as file:
features_dict = pickle.load(file)
with open('XGBoost/best_hyperparams.pickle', 'rb') as file:
hyperparams_dict = pickle.load(file)
In [62]:
%%time
july_forecasts = []
for shift in range(1,7):
data = shifted_data[shift]
print(f'Модель T+{shift}; {datetime.now().strftime("%H:%M:%S")}')
best_hyperparams[shift] = {}
for zone in zone_cluster_labels.keys():
cluster_label = zone_cluster_labels[zone]
#print(f'Зона {zone}; {datetime.now().strftime("%H:%M:%S")}')
# временной ряд региона
timeseries = data.loc[:, [zone]]
timeseries.columns = ['real_counts']
timeseries.columns.name = zone
# метаданные по региону
#zone_data = marked_data[marked_data['PU_zone_id'] == zone].shift(-shift).dropna()
# добавляем признаки
timeseries_with_features = add_regressors(timeseries, K_s, K_h, K_d)
timeseries_with_features.columns.name = zone
# нужные признаки
selected_features = features_dict[shift][cluster_label]
## подобранные гиперпараметры модели
xgb_params = hyperparams_dict[shift][cluster_label]
xgb_params.update({'max_depth': int(xgb_params['max_depth']),
'n_estimators': int(xgb_params['n_estimators'])})
#### ПРОГНОЗИРУЕМ НА ИЮЛЬ
## отделяем train и test
timeseries_train = timeseries_with_features.loc[:split_date].copy()
X_train = timeseries_train.drop(['real_counts'], axis=1)[selected_features]
y_train = timeseries_train[['real_counts']]
timeseries_test = timeseries_with_features.loc[july[0]:july[-7]].copy()
X_test = timeseries_test.drop(['real_counts'], axis=1)[selected_features]
y_test = timeseries_test[['real_counts']]
# модель
XGBR = myXGBR(**xgb_params)
XGBR.fit(X_train, y_train)
## прогноз
july_forecast = pd.DataFrame(XGBR.predict(X_test).ravel(),
index = X_test.index,
columns=['prediction'])
july_forecast = pd.concat([y_test, july_forecast], axis=1)
july_forecast['zone'] = zone
july_forecast['shift'] = shift
july_forecasts.append(july_forecast)
In [63]:
july_forecasts = pd.concat(july_forecasts, axis = 0)
In [64]:
july_forecasts.to_csv('XGBoost/july_forecasts.csv', index=True)
In [65]:
# посчитаем Q_july
Q_july = np.sum(np.abs(july_forecasts['real_counts'] - july_forecasts['prediction']))/ (R*739*6)
print('Q_july= {:.2f}'.format(Q_july))
Кажется, это не очень большая средняя ошибка прогноза в масштабах спроса на такси в Нью-Йорке.